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Abstract. We present the results of 3D Newtonian SPH 
simulations of the merger of a neutron star binary. The mi- 
croscopic properties of matter are described by the phys- 
ical equation of state of Lattimer and Swesty (LS-EOS). 
To check the model dependence of the results we vary the 
resolution (~ 21000 and ~ 50000 particles), the equation 
of state (stiff and soft polytropes), the artificial viscos- 
ity scheme, the stellar masses, we include neutrinos (free- 
streaming limit), switch off the gravitational backreaction 
force, and vary the initial stellar spins. In addition we 
test the influence of the initial configuration, i.e. spheri- 
cal stars versus corotating equilibrium configurations. The 
final matter distribution consists of a rapidly spinning cen- 
tral object with 2.5 to 3.1 M Q of baryonic mass that prob- 
ably collapses to a black hole, a thick disk of 0.1 to 0.3 M 
and an extended low density region. In the case of coro- 
tation this low density material forms spiral arms that 
expand explosively due to an increase of the adiabatic ex- 
ponent and the release of nuclear binding energy in the 
case of the LS-EOS, but remain narrow and well defined 
for the stiff polytropic equation of state. The main and 
new result is that for the realistic LS-EOS, depending on 
the initial spin, between 4 • 10~ 3 and 4 • 10 -2 M of mate- 
rial become unbound. If, as suggested, large parts of this 
matter consist of r-process nuclei, neutron star mergers 
could account for the whole observed r-process material 
in the Galaxy. 
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1. Introduction 

The question of the origin of Gamma Ray Bursts (GRBs) 
-whether they are galactic or cosmological- has been con- 
troversial during the three decades since the announce- 
ment of their detection (Klebesadel et al. 1973). In spite 
of the fact that the Burst and Transient Source Experi- 
ment (BATSE) detected bursts at a rate of around one 
per day, no clear indication of a distance scale could be 
found. It was only recently that observations of the after- 
glow of Gamma Ray Bursts in the optical (van Paradijs et 
al. 1997; Sahu et al. 1997; Galama et al. 1997; Djorgovski 
et al. 1997), X-ray (Costa et al. 1997) and radio (Frail et 
al. 1997) part of the spectrum seem to have settled their 
cosmological origin. The identification of Fe and Mg ab- 
sorption lines and the determination of the cosmological 
redshift (0.835 < z < 2.3; Metzger et al. 1997) seems to 
have definitely ruled out galactic sources of GRBs. The 
most popular model for the central engine to power cos- 
mological GRBs is the merger of two neutron stars (or 
a neutron star and a black hole; see e.g. Paczyhski 1986, 
1991, 1992; Eichler et al. 1989; Rees & Mesczaros 1992; 
Narayan et al. 1992). 

This scenario is also attractive for other reasons. By now 
five such neutron star binary systems are known (Thorsett 
1996) and no mechanism is known that could prevent the 
inspiral and final coalescence. During the last ^15 min- 
utes of its inspiral the binary system will emit gravita- 
tional waves with a frequency ranging from ~ 10 Hz to 
~ 1000 Hz. The inspiral is thus one of the prime candi- 
dates to be detected by earth bound gravitational wave 
detector facilities such as LIGO (Abramovici et al. 1992), 
VIRGO (Bradaschia et al. 1990) or GEO600 (Luck 1997). 
Another reason to study this scenario -and this will be our 
main focus here- is the question of nucleosynthesis, namely 
the production site of r-process nuclei. Despite consider- 
able efforts, it has up to now not been possible to iden- 
tify the corresponding astrophysical event unambiguously. 
Most recent nucleosynthesis studies (see Freiburghaus et 
al. 1997, 1998; Meyer & Brown 1997) raise questions con- 
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cerning the ability of high entropy neutrino wind scenar- 
ios in type II supernovae to produce r-process nuclei for 
A < 110 in correct amounts. In addition, it remains an 
open question whether the entropies required for the nu- 
clei with A > 110 can actually be attained in type II 
supernova events. Thus, an alternative or supplementary 
low entropy, low Y e r-process environment seems to be 
needed (decompression of neutron star material). 
The decompression of neutron star matter as a possible 
source of r-process nuclei was first discussed by Lattimer 
& Schramm (1974, 1976) in their study of the tidal dis- 
ruption of a neutron star by a black hole. The coalescence 
of two neutron stars has in this context been examined by 
Symbalisty & Schramm (1982) and Eichler et al. (1989). 
Detailed studies of the nuclear decompression process have 
been performed by Meyer (1989). 

Hydrodynamic simulations of coalescing neutron star bi- 
naries have been performed by various groups. The first 
calculations were done by Nakamura and Oohara (see Shi- 
bata et al. (1993) and references therein) using polytropic 
equations of state and focusing mainly on the emitted 
gravitational waves. 

Similar calculations, but using SPH, have been performed 
by Davies et al. (1994) who discussed a variety of physical 
effects related with the merging event. Zhuge et al. (1994, 
1996) focussed in their work on the gravitational wave en- 
ergy spectra dE/df. 

In a long series of papers Lai, Rasio and Shapiro exam- 
ined close binary systems. They developed an approximate 
analytical energy variational method and applied it to an- 
alyze the stability properties of binary systems and rotat- 
ing stars (Lai et al. 1993a, 1993b, 1994a, 1994b, 1994c, Lai 
1994). Two of them (Rasio and Shapiro) performed com- 
plementary SPH-simulations, where apart from stability 
questions the emission of gravitational waves was investi- 
gated (Rasio & Shapiro 1992, 1994, 1995). 
Ruffert et al. (1996, 1997) performed PPM-simulations of 
neutron star mergers, using the Lattimer/Swesty EOS and 
accounting for neutrino emission by means of an elabo- 
rate leakage scheme. They discussed gravitational wave 
and neutrino emission and made an attempt to address 
the question of mass ejection by looking at the material 
they lost from their grid. 

Rosswog et al. (1998a, b) discussed preliminary results of 
their SPH calculations concerning mass ejection and its 
implications for nucleosynthesis. 

Fully relativistic hydrodynamical calculations are begin- 
ning to yield results (see Wilson & Mathews 1995, Wil- 
son et al. 1996, Mathews & Wilson 1997 , Baumgarte et 
al. 1997). However, these are still controversial (see e.g. 
Lai 1996, Thornc 1997, Lombardi & Rasio 1997, Wise- 
man 1997, Shibata et al. 1998). 

Our focus in this paper will be on the matter that becomes 
unbound in an equal mass neutron star binary merger. In 
Sect. 2 we describe the numerical method, we discuss our 
chosen binary system parameters and the appropriateness 



of the approximations of our model. The results arc dis- 
cussed in Sect. 3. In 3.1 we describe the morphology of the 
mergers, in 3.2 the mass distribution is discussed, 3.3 deals 
with the thermodynamic properties of the merged config- 
uration, in 3.4 we discuss our neutrino treatment and 3.5 
deals with mass ejection and related questions of nucle- 
osynthesis. The results are summarized in Sect. 4. Details 
of the initial SPH-particle setup, the equation of state and 
the neutrino emission are given in the appendices. 

2. Numerical Method and Physical Parameters of 
the System 

2.1. The Numerical Method 

We perform hydrodynamic calculations using a 3D New- 
tonian SPH code based on the one used in Davies et al. 
(1994). The basic equations of our SPH formulation can 
be found in the review of Bcnz (1990), the gravitational 
forces are calculated using a hierarchical binary tree (Bcnz 
et al. 1990). Since the method is well-known and compre- 
hensive reviews exist (see e.g. Hcrnquist & Katz (1989); 
Benz (1990); Monaghan (1992); Miiller (1998)), we do not 
go into the details of SPH here and refer the reader to the 
literature. 

The Lagrangian nature of SPH is perfectly suited to the 
study of this intrinsically three dimensional problem, since 
it is not subject to spatial restrictions imposed by a com- 
putational grid. Thus, no matter can be lost from the com- 
putational domain (as for example in Eulerian methods). 
In addition no CPU resources are wasted for uninteresting 
empty regions. These zones are in other methods some- 
times filled by low density material that can possibly lead 
to artifacts. Especially, the (tiny) amounts of ejecta could 
be influenced since they have to struggle against surround- 
ing material. SPH also allows one to follow individual par- 
ticles and therefore to keep track of the history of chosen 
blobs of matter. 

The details of the initial SPH-particle setup are given in 
Appendix A. 

2.2. The Physical Parameters of the Binary System 

In the following we want to discuss the physical param- 
eters of the binary system together with the chosen ap- 
proaches and their uncertainties that will be tested. 

2.2.1. Masses 

Present state of the art nuclear equations of state allow 
gravitational masses for maximally rotating neutron stars 
of - 2.5 Mq (see e.g. Weber & Glendenning (1992)). How- 
ever, all observed neutron star masses are centered in a 
narrow band around 1.4 M© (e.g. van den Heuvel (1994)), 
a mass that might be given rather by evolutionary con- 
straints than by the stability limit of the EOS. The masses 
that are found in the five known neutron star binaries are 
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very close to each other, the maximum deviation being ~ 
4 % in the case of PSR 1913+16 (see Thorsett 1996). Con- 
sidering the observational facts we only investigate equal 
mass binary systems. We regard a baryonic mass of 1.6 
M Q , corresponding approximately to 1.4 M Q of gravita- 
tional mass, as most reasonable. The amount of mass that 
is ejected into the interstellar medium is directly related 
to the gravitational potential that has to be overcome 
by the neutron star debris. Hence we vary the neutron 
star masses, investigating also the case of 1.4 M Q of bary- 
onic matter per star. Starting with the initial separation 
ao = 45 km, which is our standard value, the inspiral will 
take longer in this case than for the heavier stars. For an 
estimate one may apply the formula for the inspiral time 
of a point mass binary (see e.g. Misner et al. 1973) 



5 c 5 



256 G 3 iiM 2 



(1) 



giving an inspiral time that is longer by a factor of ap- 
proximately (8/7) 3 ~ 1.5 for the lighter case. 

2.2.2. Gravitational Radiation Backreaction 

Since the backreaction force, resulting from the emission of 
gravitational waves, tends to circularize elliptic orbits (see 
Peters & Mathews (1963, 1964)), we start with basically 
circular orbits. Starting on the x-axis of our coordinate 
system, we give each SPH-particle a tangential velocity 
corresponding to a circular Kepler motion with a fixed 
value lo 2 — GMa 3 , where M denotes the total mass of 
the system, and add the radial velocity of a point mass 
binary in quadrupole approximation (see e.g. Shapiro & 
Teukolsky 1983) 



dr 



64 G 3 Mj 
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down) the backreaction force is switched off. 
This treatment of the backreaction force is clearly a sim- 
plification of the relevant gravitational physics. To test 
the corresponding sensitivity of our results we will inves- 
tigate also the case where this force is totally neglected. To 
push the system over the critical radius where it becomes 
dynamically unstable (for an estimate of this radius ob- 
tained by applying an energy variational method to com- 
pressible, triaxial ellipsoids obeying a polytropic equation 
of state see Lai et al. (1994a) and references therein), we 
have to begin in this case with a smaller initial separation 
(42 km). 

2.2.3. Tidal Deformation 

For the standard case of our calculation we will start with 
spherical initial configurations. However, a neutron star 
binary system that has spiraled down to a center of mass 
distance of an, = 45 km will show non-negligible tidal de- 
formations. In a crude estimate one finds for the height of 
the tidal bulge in our equal mass system h ~ (^) 3 ' R ~ 
0.5 km, i.e. parts of the surfaces of the neutron stars will 
be raised to up to half a kilometer above their spherical 
radius R. Thus starting with initially spherical stars will 
lead to oscillations. To estimate the impact of the approx- 
imation of spherical stars we compare the results with the 
case where the initial configuration has been relaxed in 
the mutual gravitational field. 

In order to obtain a hydrostatic equilibrium for the tidally 
locked, corotating binary system, we consider the coordi- 
nate system in which the neutron stars are at rest. A par- 
ticle i of mass m t is accelerated not only by the forces in 
the inertial frame F i: but also by non-inertial forces (see 
e.g. Landau 1976), so that the total force reads: 



2wxv, +wx(wxr, ) + u x r ; ,(5) 



(r = a/2, Mj is the mass of each star). As in Davies 
et al. (1994) we apply in addition to the hydrodynamic 
and gravitational accelerations the backreaction of a point 
mass binary to each SPH-particle: 
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Here 77 = ^ M a i M , E is the energy, J the angular mo- 
mentum, fj, the reduced mass and v, r, x, y, x and y 
refer to the point mass binary (not to be confused with 
the SPH-particle properties). Since this acceleration is ap- 
plied to each particle, the circulation of the fluid will be 
conserved. When the distances of the centers of mass of 
the stars from the origin are equal to one stellar radius 
(and thus the point mass approximation definitely breaks 



where us is the rotational frequency of the non-inertial sys- 
tem, i.e. the Kepler frequency of the corresponding point 
mass binary, and rj and Vi are position and velocity 
vectors of particle i in the corotating system. The first 
term in the brackets is the Coriolis-force Fcor, the second 
the centrifugal force Fcont and the last term, , results 
from the change of us due to the gravitational radiation 
backreaction force. To obtain equilibrium we apply a ve- 
locity dependent damping force, so that the total force on 
a SPH-particle i reads: 

F i t = F lygIav + F^hydro + -FY Cent + F itdl - m^v ^ (6) 

where F i grav and i*Y hydro are forces due to the presence of 
self- gravity and pressure gradients, both of which are eval- 
uated with the SPH-code. We do not have to consider the 
backreaction force here, since it is applied (in our approx- 
imation) to each SPH-particle in the same way and hence 
does not influence the shape of the stars. The Coriolis force 
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does not have to be considered here since we are inter- 
ested in the construction of hydrostatic equilibrium where 
v i = 0. As mentioned by Rasio and Shapiro (1992) the 
relaxation time r = 7 _1 should be chosen slightly over- 
critical to guarantee an efficient relaxation process (the 
typical oscillation time scale can be estimated from the 
sound crossing time t s = (Gp)~i). 

We first calculate the resulting forces on the centers of 
mass and then subtract them from the forces on individ- 
ual particles in order to remove the force on the stellar 
centers of mass: 



F 3 - F 3 

i emi 



where the center of mass forces are given by 

3 ;. 



(7) 



(8) 



Here j labels the binary stars, with Mj their masses and 
the sum loops over the particle numbers in each object. 
Because of numerical noise the centers of mass drift very 
slightly and thus are reset each dump to keep the stars 
exactly at the desired positions. 

Kepler's law lu 2 = GMa^ 3 , which we applied in the re- 
laxation process, is strictly valid only for point masses. In 
their ellipsoidal approximation for polytropic equations of 
state Lai et al. (1994a) find a modified Kepler law for close 
equal mass binaries: 



tidal 



= uj 2 (1 + 2(5), 



where 



3 2J n 
2 



'22 



'33 



(9) 



(10) 



and In = K n a 2 /5. The en are the semi-major axes of the 
ellipsoid and K n is a constant depending on the polytropic 
index n (k„ = 1 in the incompressible case, n = 0; nu- 
merical values are tabulated in Lai et al. (1993a)). We 
calculated Wtidai a posteriori using the parameters of the 
relaxed system. We found Wtidai ~ 1.015-w, which justifies 
using w instead of Wtidai m the relaxation process. 
Starting with a binary system relaxed in the way described 
above no oscillations of the stars were observed. 

2.2.4. Viscosity 

Since the neutrons and protons are likely to be in super- 
fluid states, the main contribution to the microscopic vis- 
cosity of the neutron star fluid is supposed to come from 
electron-electron scattering (Flowers & Itoh 1979). The 
simulation of the almost inviscid neutron star fluid poses 
a problem for numerical hydrodynamics since the latter is 
subject to numerical as well as to artificial viscosity. 
The SPH standard tensor of artificial viscosity (Monaghan 
& Gingold 1983; we use the standard values a = 1 and 



(3 = 2) is known to yield good results in shocks, but also 
to introduce spurious forces in pure shear flows, where a 
vanishing viscosity would be desirable. We thus decided 
to compare the results obtained using the standard vis- 
cosity with those found using a new scheme proposed by 
Morris & Monaghan (1997). The basic idea is to give each 
particle its own time dependent viscosity coefficient that 
is calculated by integrating a simple differential equation 
including a source term and a term describing the decay 
of this coefficient. The viscosity tensor then reads 



_ cci+aj 



rijV i:j < 



Ci+Cj 



hj + hj 



where = 2 a>ij, a-, 

and Hij = ^Tl' 3 ^ ■ The Ck denote the particle sound ve- 

locities, hk the smoothing lengths, and Vk positions 
and velocities, — — rj, Vij = Vi — Vj and r\ = 0.01. 
cti is calculated from 



doii 
~~dt 



ex; — a 



+ Si. 



(11) 



The first term on the right hand side causes the viscosity 
coefficient a to decay on a time scale Tj towards some 
minimum value a* that keeps the particles well ordered in 
the absence of shocks. The source term Si, leading to an 
increase of at when a shock is detected, is chosen to be 

Si = max(-(Vv)i,0). (12) 

We chose r» = ^j- with a* = e = 0.1. 

2.2.5. Spins 

The degree of synchronization of a close binary system de- 
pends decisively on the viscous dissipation rate (compared 
to the inspiral time). Bildstcn & Cutler (1992), Kochanck 
(1992) and Lai (1994) found any reasonable assumption 
on the microscopic viscosity of the neutron star matter to 
be orders of magnitude too low to lead to a tidal locking of 
a binary system within the inspiral time. Hence, it seems 
very unlikely that a corotating system is the generic case. 
We suspect the viscosity to be also too low to lead neces- 
sarily to an alignment of the stellar spins with the orbital 
angular momentum. However, for reasons of simplicity we 
will only treat aligned systems in this work. 
For purely numerical reasons, i.e. to minimize the effects 
of viscosity during the inspiral and the shear motion at 
the contact surface of both stars during the merger, we 
focus here mainly on synchronized, "corotating" systems. 
This might also be useful as some kind of benchmark for 
future simulations with a more sophisticated incorpora- 
tion of general relativistic effects (for a general relativistic 
simulation of a corotating inspiral see e.g. Baumgarte et 
al. 1997). 

However, to explore the dependence of our results on the 
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initial neutron star spins, we also investigate the cases 
where the stars have no initial spin (in a space- fixed frame) 
and the case where they are spinning against the orbital 
motion with a period equal to the orbital period. Start- 
ing with 45 km of initial separation, this corresponds to 
a spin period of 2.91 ms. We expect all mergers to lie in 
this range of initial spins. 

2.2.6. Equation of state 

The equation of state is a major uncertainty in all neu- 
tron star related calculations. Our code is coupled to the 
physical equation of state (EOS) for hot and dense mat- 
ter of Lattimer and Swesty (1991). This EOS models the 
hadronic matter as consisting of an average heavy nucleus 
(representing all heavy nuclei), alpha particles (represen- 
tative for an ensemble of light nuclei) and nucleons outside 
nuclei. Since further hadronic degrees of freedom, such as 
hypcrons, pions, kaons and the quark- hadron phase tran- 
sition are disregarded, we use the lowest available nuclear 
compressibility, K — 180 MeV, to mimic the possible 
softening of the EOS due to the appearance of "exotic" 
matter at higher densities. Since the chemical composi- 
tion of matter is calculated for the thermodynamically 
most favourable state, effects of recombinations of nucle- 
ons into nuclei are taken into account automatically. For 
further details concerning the LS-EOS see Appendix B. 
To test the dependence of our results on the underlying 
EOS we also use a polytropic EOS where the exponent V 
was fitted to the central part of our 1.6 M Q neutron star 
obtained with the LS-EOS. We obtained a good fit with 
T = 2.6. The pressure p can be written in terms of density 
p and specific internal energy u p 



p=(T-l)pu 



v 



(13) 



We assign the specific internal energy of our polytrope, 
u p , in a way that the pressure distribution of the LS-EOS 
star is reproduced: 



u p (r) 



PLsjr) 

(r-i)p(rf 



(14) 



where PL.s( r ) an d p(r) are the values of pressure and den- 
sity at radius r in the star. Using these profiles a spherical 
equilibrium configuration was constructed following the 
steps described in Appendix A. 

To test the sensitivity on the adiabatic exponent we also 
considered a test case of a much softer EOS (F = 2.0) 
where the polytropic constant K (p — Kp T ) was adjusted 
in a way that the radius of the LS-EOS star was repro- 
duced. 

2.2.7. Temperatures 

When the neutron stars of a binary system have reached 
a separation of a few stellar radii (typical time scales are 



of the order of 10 8 years) they should be essentially cold 
(T < 10 6 K). When the stars come closer and tidal inter- 
actions become important, the stars might be heated up 
again by viscous processes acting in the tidal flow. How- 
ever, these temperatures are estimated to be only of the 
order of 10 8 K (Lai 1994). 

With respect to the thermodynamic properties, the accu- 
rate determination of low temperatures in the high den- 
sity regime (say above 10 14 gcm -3 ) is very difficult. Since 
the temperature is a very steeply rising function of the 
specific internal energy u (our independent variables are 
u, the electron fraction Y e and the density p), numerical 
noise in u in the degenerate regime may lead to substan- 
tial temperature fluctuations. However, as long as there 
are no additional physical processes involved that depend 
on temperature, this does not influence the dynamical evo- 
lution of matter. 



2.2.8. Neutrino emission 

The merger of two neutron stars will be accompanied by 
strong neutrino emission (see e.g. Ruffert et al. (1996, 
1997)). To set a limit on the maximum influence of the 
cooling by neutrino emission on the ejection of mass, we 
also explore the extreme case where we assume that all lo- 
cally produced neutrinos stream out freely without any in- 
teraction with overlying material (which also implies that 
a possible quenching of neutrino rates due to neutrino 
final-state blocking is ignored), thereby reducing the local 
thermal energy content. This scheme is clearly crude and 
will largely overestimate the emission from the hot, dense 
central regions. However, we are here only interested in 
the possible influence on the amount of ejected material 
and the free streaming limit will give us a secure upper 
limit of the effects resulting from neutrino emission. 

Looking at typical densities and temperatures encoun- 
tered in our scenario, one can conclude that the only pro- 
cesses giving possibly non-negligible contributions to the 
neutrino emission in addition to the electron and positron 
capture on nucleons, are the plasma and the pair processes 
(for a comparison of the importance of different neutrino 
processes as a function of p, Y e and T see e.g. Schinder et 
al. (1987), Itoh et al. (1989) and Haft et al. (1994)). 
We thus include the following processes as sources of en- 
ergy loss: 



e + p 
e + + n 
e+ + e 

1L,T 



n + v e 
p +P e 
Vi + v l 



(15) 
(16) 
(17) 
(18) 



where the last process is the decay of a plasmon, i.e. the 
decay of both longitudinal and transverse electromagnetic 
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excitations of the medium. Due to the very high tempera- 
tures, we ignore the electron mass m e and the mass differ- 
ence Q of neutron and proton for our treatment of electron 
and positron captures. The energy emission rates due to 
these processes are discussed in Appendix C. For the en- 
ergy emission rates of pair and plasma neutrino processes 
we use the fit formulae of ftoh et al. (1989) and Haft et 
al. (1994). The neutrino emission then appears as a sink 
in the energy equation of each particle: 

(hH _ fdUi\ QEC + QPC + Qpair + Qplas ,-.„x 

dt { ill ).. Pi 

where (^ i ) contains the usual terms from pdV -work and 
viscous heating (see e.g. Benz (1990)). The enormous de- 
pendence of the rates on temperature (QecQpc ~ T 6 ; 
Qpair,Qpias ~ T 9 ) gave rise to numerical instabilities in 
the degenerate regime, where temperature fluctuations 
can be appreciable due to noise in the specific internal 
energy. To overcome this problem the emission rates of 
particle i were calculated from mean temperatures T; t = 
N^ 1 J2f l Tj , where Ni denotes the number of neighbours 
of particle i, rather than from the individual particle tem- 
peratures. 

3. Results 

After having discussed in Sect. 2 the numerical method 
(including the initial setup) and all possible sensitivities 
and uncertainties, we have performed in total 11 numeri- 
cal simulations (run A to K), to test all of these aspects. 
In this way we intend to analyze the robustness of the 
results from the (corotating) standard case and to un- 
derstand what variation of these results can be expected 
for different system parameters and physics input. The 
specifics which characterize each of these runs are given 
in Table 1. The following discussion of the results will re- 
fer to these table entries. In all of these runs we follow the 
evolution of matter for 12.9 ms. 

3.1. Morphology 

Figs. 1 to 4 show the particle positions of runs A, C, I and 
J. The different colors in Fig. 1 label those particles that 
end up (last dump at t — 12.9 ms) in the central object 
(black), the disk (green), the tails (blue) and those that 
are unbound at the end of the simulation (red). Run A 
stands representative for all the corotating runs using the 
Lattimer/Swesty EOS (the particle plots of runs A, B, D 
and E are practically indistinguishable; runs F, G and H 
differ only slightly). The plots only show the first half of 
the simulation. In the second half practically no additional 
angular motion is visible, matter just expands radially up 
to ~ 900 km for all runs apart from the one without initial 
stellar spins. In this case the outermost layers only extend 
out to <~ 600 km from the center of mass. 



In all corotating runs those parts having the largest dis- 
tance from the common center of mass are driven into 
thick spiral arms. Those get "wrapped up" around the 
central object to form a ring-like structure (see Figs. 1 and 
2). In the further evolution this ring contracts to form a 
thick disk, while the matter in the spiral arms continues 
to wind up around the coalesced, massive object in the 
center of mass. 

The corotating runs show, apart from the "main" spiral 
arms, also a very narrow "side" spiral structure (see, for 
example, the second panel in Fig. 2; in outlines, this struc- 
ture is also visible in run I), that emerges directly after the 
stars have come into contact. The side spiral arm matter 
comes from the "contact zone" (see panel one in Fig. 1), is 
then strongly compressed and finally, in a kind of "tube- 
of-toothpaste" -effect ejected into two narrow jets in oppo- 
site directions. While matter in the main arms is smoothly 
decompressed due to gravitational torques, matter in the 
side spiral arms gets compressed and heated up. The pres- 
sure responsible for the formation of the side spiral arms 
is of thermal nature. This is suggested since no such spiral 
arms emerge in run G, where neutrino emission leads to a 
very efficient cooling. In an additional test run, where the 
viscous pressure contribution to the energy equation was 
disregarded, also no side spiral arms occurred. 
In the end three well-separated structures are visible: an 
oblate, massive central object, a thick surrounding disk, 
and long spiral arms. Going to lower total angular mo- 
mentum, i.e. to runs I and J, the differences between these 
structures get washed out. 

It is interesting to note the dependence of the spiral arms 
on the equation of state: while they are thick and pumped 
up in all corotating runs using the LS-EOS, they remain 
narrow and well-defined in the run with the stiff polytropic 
EOS (r = 2.6, run C; somewhat thicker for Y = 2.0, run 
K). This expansion is a result of an increase of the adi- 
abatic index and the energy release due to the recombi- 
nation of the nucleons into heavy nuclei. Fig. 5 shows the 
total amount of nuclear binding energy present according 
to the LS-EOS: 



S nuc ,tot = J2 mi ^ B h + Y a B a ], (20) 



where the summation ranges over all SPH-particles, m t 
denotes the particle mass, Yh and Y a the abundances of 
heavy nuclei and alpha particles, and Bh and B a the nu- 
clear binding energies. Around 10 50 erg of nuclear binding 
energy are deposited (i.e. present in the last dump) in the 
spiral arms containing only ~ 0.07 M Q (see Table 2). 
The degree of conservation of energy and angular momen- 
tum can be inferred from Fig. 6. Both are conserved to 
approximately 3 • 10~ 3 . 
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Table 1. Summary of the different runs: 1: corotation; 2: no initial spins; 3: spin against the orbit; LS: Lattimer and Swesty 
(1991), P: polytropic EOS ; MG: Monaghan and Gingold (1983); MM: Morris and Monaghan (1997) 



run 


spin 


EOS 




viscosity 


backreaction 


tidal def. 


V 


M [Mol 


$ part. 


do [km] 


A 


1 


LS 




MG 


yes 


no 


no 


1.6 


20852 


45 


B 


1 


LS 




MG 


yes 


yes 


no 


1.6 


20852 


45 


c 


1 


p, r = 


2.6 


MG 


yes 


no 


no 


1.6 


20852 


45 


D 


1 


LS 




MM 


yes 


no 


no 


1.6 


20852 


45 


E 


1 


LS 




MG 


yes 


no 


no 


1.6 


49974 


45 


F 


1 


LS 




MG 


yes 


no 


no 


1.4 


20852 


45 


G 


1 


LS 




MG 


yes 


no 


yes 


1.6 


20852 


45 


H 


1 


LS 




MG 


no 


no 


no 


1.6 


20852 


42 


I 


2 


LS 




MM 


yes 


no 


no 


1.6 


49974 


45 


J 


3 


LS 




MM 


yes 


no 


no 


1.6 


49974 


45 


K 


1 


p,r = 


2.0 


MG 


yes 


no 


no 


1.6 


20974 


45 



Table 2. Masses of the different morphological regions, M . o . is the minimum gravitational mass (see Eq. (22)), a is the 
relativistic stability parameter (see text). 



run 


M c .o. [M ] 


Mc.o. [M ] 


M disk [M ] 


Mfii [M©] 


a 


A 


2.81 


2.32 


0.32 


0.07 


0.54 


B 


2.81 


2.32 


0.32 


0.07 


0.54 


C 


2.78 


2.29 


0.35 


0.08 


0.51 


D 


2.81 


2.32 


0.32 


0.07 


0.51 


E 


2.81 


2.32 


0.32 


0.07 


0.54 


F 


2.46 


2.07 


0.27 


0.07 


0.56 


G 


2.88 


2.36 


0.25 


0.07 


0.58 


H 


2.75 


2.27 


0.32 


0.12 


0.49 


I 


3.05 


2.47 


0.11 


0.04 


0.62 


J 


3.14 


2.54 


0.06 


8 ■ 10~ 3 


0.61 


K 


2.76 


2.28 


0.40 


0.04 


0.53 



3.2. Mass distribution 

Density contours (cuts through the y-z-plane at the last 
dump, t=12.9 ms) of run E (representative for the corotat- 
ing runs), run G (to see the influence of the cooling by neu- 
trino emission) , run I and run J (to see the influence of the 
initial spins) are shown in Figs. 7 and 8. The distribution 
of mass with density and radius is shown in Figs. 9 and 10. 
The sharp bends in Fig. 9 indicate the separations of the 
different morphological regions. The amounts of mass in 
the central object, the disk and the tails (collectively used 
for the low-density material outside the disk) can be in- 
ferred from this plot, the corresponding masses are listed 
in Table 2. The runs A, B, D and E practically coincide 
in Fig. 9, suggesting that the mass distribution is insen- 
sitive to our approximation of initially spherical stars, to 
the details of artificial viscosity, and also shows that the 
resolution in our standard run is sufficient to describe the 
mass distribution properly. The temporal evolution of the 
masses and densities in the three regions (run A) can be 
inferred from Fig. 11 (Fig. 9 just corresponds to the last 



temporal slice of Fig. 11). In the beginning practically all 
the mass (= 3.2 M ) has a density above 10 14 gcm -3 (the 
visible oscillations result from the initially spherical stars). 
Around t= 2.5 ms a strong expansion sets in (less mass 
has a density above, say, 10 gem -3 ). Then, around t= 4.5 
ms, the three regions (central object, disk and tails) be- 
come visible. At around t= 5.5 ms, when the wrapped-up 
spiral arms are contracted by gravity, a recompactifica- 
tion corresponding, to the "bump" in Fig. 11 sets in (to 
localize it in time and density the contour lines of 2.8, 2.9, 
3.0 and 3.1 M Q are projected on the log(/?)-t-planc). This 
recompactification is accompanied by a reheating (this is 
also reflected in the neutrino emission of the disk, where a 
strong increase is visible, compare panel 2 in Fig. 22) that 
finally reexpands the disk. When our simulation stops, the 
disk has a density range from ~ 10 10 — 10 12 gem -3 (com- 
pare also to the contour plots 7 and 8), densities above 
this range are associated with the central object, densities 
below with the tails. 

The overall mass distribution is shown in Figs. 12 and 13. 



8 



S. Rosswog et al.: Mass ejection in neutron star mergers 



3.2.1. The central object 

Figs. 14 and 15 show density contours of runs E, C, I 
and J. The central objects (especially for runs A to H) 
consist of a core with densities above 10 14 gcm~ 3 and a 
diffuse edge ranging from 10 14 to ~ 10 125 gcm~ 3 forming 
a kind of a "hot skin" around the core (see the smooth 
transition from the central object to the disk in Fig. 9). 
In the run including neutrino energy losses (run G) there 
is a very sharp transition corresponding to a huge density 
gradient from the central object to the disk, suggesting 
that the edge diffuscncss of the central object in the other 
runs emerges from thermal pressure (also visible in the 
temperature, Fig. 19). Because of the absence of this pres- 
sure, matter of the order 0.07 M Q , otherwise located in the 
disk, has fallen onto the central object increasing its mass 
to 2.88 M Q . (This hot skin might be an artifact of the 
viscosity scheme, since it less pronounced in run D). The 
masses of the final central objects range from ~ 2.5 M Q to 
~ 3.1 M Q (see Table 2) and are thus well above the maxi- 
mum precisely known neutron star (gravitational) mass of 
1.44 M for the pulsar of the binary system PSR 1913+16 
(Taylor 1994; there are, however, objects, whose mass er- 
ror bars reach up to 2.5 M Q ; see Prakash (1997b) and 
references therein) . Theoretical calculations for maximally 
rotating neutron stars, using a large set of 17 (11 relativis- 
tic nuclear field theoretical and 6 non-relativistic potential 
models for the nucleon-nucleon-force) nuclear equations of 
state (Weber & Glcndcnning 1992) find maximum values 
for the gravitational masses of ~ 2.5 M Q . Since these val- 
ues refer to gravitational rather than to baryonic masses 
(which we are referring to due to the Newtonian character 
of our calculations) we have to estimate the gravitational 
masses of our central objects. The relation obtained by 
Lattimer and Yahil (1989) for reasonable uncertainties in 
the equation of state reads 



1 -"Hal 



10 53 erg. 



(21) 



This gives the minimal gravitational mass M 9 of the bary- 
onic mass M b as 



5.420. [-1 + 



1 + 0.369 



M b . 



(22) 



which is typically 0.4 - 0.5 M Q smaller than the baryonic 
mass. The corresponding values for our central objects are 
shown as entry 3 in Table 2). These values are still very 
high, but there are equations of state for which fast rotat- 
ing neutron stars in this mass range are stable. Thermal 
pressure, which is disregarded in the work of Weber and 
Glcndenning since they assume old neutron stars, could 
further stabilize the central object on a cooling time scale. 
The above quoted masses refer to neutron stars rotat- 
ing with the maximum possible velocity. To estimate the 
influence of rotation for our calculations we plot in Fig. 



17 the tangential velocities in the central objects of run 
E, I and J and compare them to the Kepler velocity 
vk = {GM^/r) 1 / 2 , i.e. the velocity where the centrifugal 
forces balance the gravitational forces (M(r) is the mass 
enclosed in the cylindrical radius r). In all cases the veloc- 
ities are clearly below the Kepler velocity and we thus do 
not suppose rotation to play an important role for stabi- 
lization (in this Newtonian consideration) . In their general 
relativistic analysis Weber and Glcndcnning find the max- 
imum possible rotation frequency fix ~ 1-2 • 10 4 s _1 for 
the equations of state referred to above corresponding to 
a maximum stable mass of M « 2.5 M (Weber & Weigel 
1989; Weber et al. 1991). This is larger than our rotation 
frequencies by a factor of about two. 
In the general relativistic case the stability support from 
rotation is determined by the dimensionless parameter 



Jc 
GM 



2- (Stark & Piran 1985). We give a for our central 
objects in column six in Table 2. We find values of ~ 0.5 in 
agreement with previous simulations (Shibata et al. 1992, 
Rasio & Shapiro 1992, Ruffert et al. 1996). This is well be- 
low the critical value a cr u ~ 1, meaning that the central 
object cannot be stabilized against collapse by rotation. 
Prakash et al. (1995) studied the influence of the compo- 
sition on the maximum neutron star mass. Their general 
result is that trapped neutrinos together with nonleptonic 
negative charges (such as T,~ hyperons, or d and s quarks) 
lead to an increase of the maximum possible mass (in con- 
trast to nucleons-only matter). Thus, the collapse of a star 
with almost the maximum mass could be delayed on a 
neutrino diffusion time scale. As an estimate of this time 
scale for our central objects, we use R « 15 km (center to 
pole distance, see Fig. 7) and a typical neutrino energy of 



F 5 (0) 
F 4 (0) 



(E Va >«T 



F 5 (0) 
F 4 (0) 



50 MeV, where T« 10 MeV and 



5 has been used. This gives a mean free path of 



A w — « 0.75 m, 
no~ 



(23) 



where the baryon number density n corresponding to 5 • 
10 14 gcm~ 3 and a neutrino nucleon scattering cross section 

(j = ±er (^rrs) with °o = 1.76- 10~ 44 cm 2 (see Tubbs & 

41 \ Ttl G C J 

Schramm 1975) has been used. An estimate for the delay 
time scale is then given by 



3R 2 
Ac 



3 s. 



(24) 



Thus if there were nonleptonic negative charges present in 
the central object, the collapse could be delayed for a few 
seconds. 

To summarize the stability question of the central object: 
if all the stabilizing effects mentioned above (EOS, rota- 
tion, thermal pressure and exotic matter) should conspire, 
the central object could (at least in some cases) be sta- 
bilized against gravitational collapse. However, we regard 
this possibility to be fairly unlikely. 
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If the central object collapses to a black hole, the question 
arises of how much mass has enough angular momentum 
to avoid being swallowed. For a simple estimate we as- 
sume that the specific angular momentum of a particle 
must be larger than the one of a test particle with Ke- 
pler velocity at the marginally stable circular direct or- 
bit of a Kerr black hole (see Bardeen et al. 1972). Thus, 
k = \r, x Vi\ > R ms VK{Rms) = (GM c . ..R m s) _1/2 , where 



i? ms = ^M c . . {3 + z 2 - [(3 - «i)(3 + Z! + 2z 2 )] 1 /2 }(25) 
with 



zi= 1+11- 



,2 \ 1/3 



M 2 

A v± C . O . 



1 + 4^-1 + 



Af c . . 



1/3 



1/3' 



and 



■Al 



c.o. , 2 



M 2 



1/2 



(26) 



(27) 



Z c .o. denotes the specific angular momentum of the central 
object. With this estimate we find that the central black 
hole will be surrounded by a disk of ~ 0.3 M Q for initial 
corotation (~ 0.1 and ~ 0.06 for the cases of no initial 
stellar spins and spins against the orbit). 

3.2.2. The disk 

In all runs the central object is finally surrounded by a 
thick disk (see Figs. 7 and 8). In the case of initial corota- 
tion the disks contain around 0.3 M Q (0.27 M in run F; 
0.25 in run G, where 0.07 M have fallen onto the central 
object; in run H, neglecting the backreaction force leads 
to an impact with a higher angular momentum and thus 
more mass is expelled into the tails; the disk, however, also 
contains 0.3 M Q ). In all cases, apart from run J, empty 
funnels form above the poles of the central object. These 
are strongly enlarged in run G, where the efficient cooling 
by neutrinos leads to strongly reduced thermal pressure 
contributions, i.e. to a softening of the EOS that makes 
matter more prone to the centrifugal forces thereby lead- 
ing to a flattening of the central object as well as of the 
disk. These funnels have the appealing feature that here 
large gradients of radiation pressure can be built up which 
can lead to two well-collimated jets in opposite directions, 
pointing away from the poles. These funnels would be an 
ideal site for a fireball scenario (Goodman 1986; Shemi & 
Piran 1990; Paczynski 1990; Piran & Shemi 1993) since 
they are practically free of baryons. A baryonic load as 
small as 10 -5 M Q could prevent the formation of a GRB. 
Our present resolution (the lightest SPH-particlc has a 
mass of a few times 10~ 5 M ) is presently still too low to 
draw conclusions on this point. 



3.2.3. The tails 

Despite considerable morphological differences -exploding 
tails with the Lattimer/Swesty EOS, well-defined, narrow 
tails for the polytrope- the amount of mass in the tails is 
rather insensitive to the EOS (see Table 2). It is mainly 
dependent on the total amount of angular momentum dur- 
ing the impact, thus leading to the most massive tails in 
run H (no angular momentum lost in gravitational waves) 
with a decreasing tendency going from the standard run 
to runs I and J. 

3.3. Temperatures and vortex structures 

Figure 18 shows the maximum temperatures during the 
different runs. The curves in the upper panel refer to the 
maximum temperature of a single particle, the ones below 
to the maximum smoothed temperature (T)j, 



(T)i 



„• Pi 



3b 



tlij ) , 



(28) 



where Tj are the temperatures, rrij the masses, pj the den- 
sities, W the spherical spline kernel (see e.g. Benz 1990) 
and hij the arithmetic mean of the smoothing lengths of 
particles i and j. 

The smoothed (the particle) temperatures reach peak val- 
ues of about 50 (80) MeV for run J where the most violent 
shear motion is present, around 45 (70) MeV in run I and 
about 30 (50) MeV in the corotating runs. The resolution 
and viscosity seem to be of minor importance for the tem- 
perature calculation. Starting with spherical stars leads to 
an increased temperature since oscillation energy is trans- 
formed into heat (see also Table 3). In the corotating runs 
the hot band that forms at the contact surface dissolves 
into two hot spots (see Fig. 19). This structure develops 
further to finally form a hot, s-shaped band through the 
central object. For an understanding of these structures 
we plotted in the right columns of Figs. 19 to 21 the pro- 
jections of those particles that are contained within a thin 
slice (|,z| < 0.5 km). The projected particle positions of 
star one and two are marked with different symbols. One 
sees the formation of two macroscopic vortices that can 
be identified with the hot spots. The panels in the sec- 
ond line of Fig. 19 show patterns that are typical for 
Kclvin-Hclmholtz instabilities (see e.g. Drazin and Reid 
1981). The basic properties of this process, the formation 
of a hot band in the contact region, separation into two 
(hot) vortices and the final s-shaped hot band, are un- 
affected by resolution and the change from the standard 
SPH-viscosity to the new scheme. The finger like struc- 
tures that extend into the matter of each star (last panel 
Fig. 19) are just broader with lower resolution (run A). 
With the new viscosity scheme substructures form along 
the "fingers" leading to a more fractal appearance of the 
line separating the material of both stars. This is a result 
of a lower viscosity which leads to a dissipation of the en- 
ergy of turbulent eddies on smaller scales A„ « LR~ 4 / 3 
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Table 3. Kinetic and thermal energies (in erg) in the different morphological regions. 



run -firkin, c.o. ^kin, disk -^kin, tail -^thcrm, c.o. £^thcrm, disk -^thcrm, tail 



A 


1.0 


10 53 


1.6 


10 52 


1.4 


10 51 


2.3 


10 52 


4.6 


10 51 


3.9 


10 49 


B 


9.9 


10 52 


1.5 


10 52 


1.5 


10 51 


2.2 


10 52 


4.1 


10 51 


3.3 


10 49 


C 


8.3 


1Q 52 


1.6 


10 52 


1.1 


10 51 


1.3 


10 52 


1.9 


10 51 


4.5 


10 49 


D 


1.0 


10 53 


1.7 


10 52 


1.5 


10 51 


2.3 


10 52 


4.6 


10 51 


3.8 


10 49 


E 


1.0 


10 53 


1.8 


10 52 


1.5 


10 51 


2.3 


10 52 


4.4 


10 51 


4.0 


10 49 


F 


7.6 


10 52 


1.2 


10 52 


1.5 


10 51 


1.8 


10 52 


3.5 


10 51 


3.9 


10 49 


G 


1.2 


10 53 


1.4 


10 52 


1.4 


10 51 


1.1 


10 52 


1.9 


10 51 


3.5 


10 49 


H 


8.5 


10 52 


1.6 


10 52 


2.7 


10 51 


2.1 


10 52 


4.8 


10 51 


2.2 


10 50 


I 


1.3 


10 53 


5.1 


10 51 


3.4 


10 50 


4.2 


10 52 


1.9 


10 51 


1.3 


10 49 


J 


1.4 


10 53 


1.9 


10 51 


2.1 


10 50 


7.7 


10 52 


1.1 


10 61 


2.7 


10 49 


K 


9.1 


10 52 


2.0 


10 52 


4.7 


10 50 


2.1 


10 53 


3.5 


10 51 


5.0 


10 49 



(see e.g. Padmanabhan (1996)), where L is is a macro- 
scopic scale and R the Reynolds number. 

In run I three macroscopic vortices, a large central one 
and two smaller ones to the left and right, form along the 
contact surface. The two smaller vortices get attracted 
by the central one and fuse on a time scale of approxi- 
mately one millisecond. Approximately eight milliseconds 
after impact the material of both stars is well mixed and 
wrapped up around the central vortex. 
In the case with maximum shear motion, run J, we find 
one large vortex in the centre and several smaller ones 
along the contact surface. The smaller ones on each side 
of the center merge so that there are in total three such 
macroscopic vortices. The outer vortices move towards the 
origin and finally merge with the central one. The mate- 
rial of both stars gets mixed turbulcntly on a time scale 
of a few milliseconds. 

Ruffcrt et al. (1996) also find two macroscopic vortices for 
the corotating case. However, they find in their simula- 
tion that the vortices dissolve practically independently 
of the initial spin state into a ringlike structure, while our 
hot spots develop into an s-like shaped hot band along 
the line separating the matter of the different stars. Since 
Ruffcrt et al. do not find a particular influence of the shear 
motion on the growth time scale and since their results for 
different initial spins look very similar, they question on 
the Kelvin-Helmholtz picture and propose an alternative, 
macroscopic explanation for the flow pattern. A quantita- 
tive analysis of the incompressible, inviscid case in terms 
of linear normal mode analysis yields a growth rate for an 
unstable mode of wavelength A (for the case p\ = p2', see 
e.g. Padmanabhan (1996)) 

a^Imiuj)- 1 = —. (29) 

This implies that the shortest wavelengths will grow 
fastest and that the perturbation should grow faster with 
larger shear velocity. 



The shortest wavelengths to grow are determined by our 
numerical resolution. Let us assume that this length scale 
corresponds to the typical distance over which neighbours 
can interact, i.e. 2hi, where the smoothing length hi is to 
be taken in the shear region. Then we find A ~ 2hi ~ 5 
km. We then look at the shear velocities v. Here we find 
by looking at relative velocity projections along the shear 
interface values of ~ 0.1c for the corotation, ~ 0.3c for run 
I and ~ 0.5c for run J. Thus typical growth time scales in 
the different runs should be a a ~ 5 • 10~ 5 s, 07 ~ 2 • 10 -5 
s and a j ~ 1 ■ 10~ 5 s, where the subscripts label the runs 
(A is representative for corotation). This means that the 
perturbations have enough time to grow into the macro- 
scopic regime on a dynamical time scale of the system 
(milliseconds) . The growth times scale approximately like 
1:2:5 and this is approximately what is seen in our sim- 
ulations. Since in our calculations a strong dependence 
of the growth time scale on the shear motion, consistent 
with the Kelvin-Helmholtz time scales, is encountered and 
the vortex structures for different initial spins are clearly 
different, we do not see the necessity of an alternative ex- 
planation. 

3.4- Neutrino emission 

The main reason for an additional run with a very sim- 
ple inclusion of neutrinos is to investigate whether neu- 
trino emission has a noticeable effect on the amount of 
ejected mass. To estimate the appropriateness of our neu- 
trino treatment for the different regions we estimate typi- 
cal neutrino diffusion time scales. A typical neutrino diffu- 
sion time scale for the central object is of the order of a few 
seconds (see discussion above). Thus our neutrino treat- 
ment is definitely inappropriate there (typical time scales 
are milliseconds). Applying the same formula (24) for the 
disk, using E v = 10 McV, p = 10 n gcnr 3 and R = 60 
km, we find r^dmk ~ 4 • 10~ 4 s, which is comparable to the 
dynamical time scales. Thus at least in the outer regions 
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of the disk, the free streaming approximation might be 
justified. In the tails free streaming neutrinos are clearly 
a good approximation (apart from, perhaps, the very first 
moments after impact). In Fig. 22 we compare the amount 
of nuclear binding energy present (see Eq. (20)) to the 
amount of energy radiated in neutrinos 

E v (t) = f {L EC (t')+L PC {t')+L pail (t')+L p ^(t'))dt',(30) 
Jo 

where L\ denotes the neutrino luminosity of process A. 
The L\ are given by 

= ? ~~piW mh (31) 

where Q\a is the energy emission rate of process A and 
particle i (in erg s _1 cm -3 ). Until the end of the simulation 
the energy lost in neutrinos exceeds that gained from nu- 
clear processes by two orders of magnitude. This is mainly 
due to the central object, where unphysical amounts of 
neutrinos are emitted. In reality we expect the disk to be 
the dominant source of neutrino emission. When the disk 
reaches the above mentioned reheating phase the neutrino 
emission exceeds the nuclear energy by about one order of 
magnitude (see panel two in Fig. 22.) 
The most important result is that in the spiral arms more 
energy is gained by nuclear processes than is lost by neu- 
trino emission (in spite of our overestimate). Such conclu- 
sions were already reached in Davies et al. (1994) for low 
density regions. This is a crucial point since it shows that 
our results concerning mass ejection are independent of 
the neutrino treatment. 

3.5. Ejected mass and nucleosynthesis 

We regard a particle to be unbound if the sum of its en- 
ergies -macroscopic as well as microscopic- is positive, i.e. 
if 

Epot + E kin + E int > 0. (32) 

For E int we count all kinds of internal energies apart from 
the (negative) nuclear binding energies, i.e. 

-Eint = E 1 + E e + En, (33) 

where the indices 7, e and N denote photons, electrons 
and nucleons. All these terms are taken from the Lat- 
timer/Swesty EOS (for the polytropic case we just use 
-Eint = Uirrii). We do not consider nuclear binding ener- 
gies, since we regard an isolated nucleus at T = (the 
only internal energy comes from nuclear binding) with 
Epot + Ekin > as unbound. However, near the end of 
the evolution E- lnt is negligible and does practically not 
influence the total amount of ejected material. 
This criterion can be cross-checked by adopting a simple 
model, where we regard the particles as free point masses, 
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i.e. we neglect hydrodynamic forces resulting from pres- 
sure gradients and disregard internal degrees of freedom 
(internal energies). The particles are supposed to move 
on Kepler orbits around a point mass in the origin with 
M = M c . . + Mdisk (M c . . is the mass of the central ob- 
ject, see Table 2), which is large compared to the particle 
masses r/ij, ^ < 1- Under these assumptions the numer- 
ical eccentricities of the orbits are given by 

= f+Wp- (34) 

where Ei is the sum of the particle's kinetic and potential 
energy and Ji its angular momentum. We generally find a 
good agreement of both criteria, with deviations lying in 
the range of 1%. 



Table 4. Amount of mass that is unbound at the end of the 
simulation 

run remark m osc [M Q ] ( # part.) 



A 


'standard' (see text) 


3.1 


10" 


-2 


(295) 


B 


tidally deformed stars 


3.4 


10" 


-2 


(324) 


C 


polytropic EOS, T = 2.6 


1.5 


10" 


-2 


(131) 


D 


lower viscosity 


3.3 


10" 


-2 


(318) 


E 


double particle number 


3.2 


10" 


-2 


(844) 


F 


1.4 M Q 


3.6 


10" 


-2 


(373) 


G 


neutrinos 


3.1 


10" 


-2 


(291) 


H 


no backreaction 


5.4 


10" 


-2 


(523) 


I 


no spins 


8.7 


10" 


-3 


(208) 


J 


spins against orbit 


5.2 


10" 


-3 


(130) 


K 


polytropic EOS, V = 2.0 




0(0) 





All corotating runs using the LS-EOS eject around 
3- 1CU 2 M©, about twice as much as the run using the stiff 
polytropc. This is caused by a variation of the adiabatic 
exponent and the formation of nuclei when matter is de- 
compressed and thereby releases the gained nuclear bind- 
ing energy. For reasons of illustration we plot in Fig. 23 the 
amount of unbound mass versus the time and compare this 
with the amount of nuclear binding energy present in the 
mass that ultimately escapes (see Fig. 24) . The flatness of 
the curves in Fig. 23 indicates that not much more mass 
will be ejected during the further evolution. The rise in 
the adiabatic exponent and the deposition of a few times 
10 49 erg in ~ 3 • 1CT 2 M Q (see Table 4) leads to an explo- 
sive expansion of the spiral arm tips, thereby supporting 
the ejection of mass. Tidal deformation (run B) leads to 
an increase of the ejected mass since the system contains 
more angular momentum due to its elongated shape. As 
expected, the 1.4 M© run ejects more mass since: the grav- 
itational potential to be overcome is shallower than in the 
1.6 M© case. We suspect this to be true also in the general 
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relativistic case where the less massive stars have larger 
radii and their outer parts therefore contain more angular 
momentum. The realistic reduction of viscosity also tends 
to increase the amount of ejecta. The inclusion of neutri- 
nos does not alter the results concerning unbound mass 
(the emitted energy in neutrinos is approximately one or- 
der of magnitude lower than the released nuclear binding 
energy for the ejected matter, see panel three in Fig. 22). 
The basic intention of run H was to test the sensitivity of 
the amount of ejected mass on the details of the treatment 
of the gravitational radiation backreaction force. Here al- 
most twice as much matter is ejected (since no angular 
momentum is lost in gravitational waves), indicating that 
our results could be influenced by our simplified treatment 
of this force (see Eqs. (3) and (4)). In the runs with lower 
initial angular momentum (run I and J) the amount of 
ejecta is substantially lower, indicating a very strong de- 
pendence on the initial spins. The sensitivity to the EOS 
is underlined by the fact that in the run with the soft poly- 
trope (r = 2.0) no resolvable amount of mass is ejected, 
in agreement with the result of Rasio and Shapiro (1992). 
The amount of r-process material that could be formed in 
this merging scenario is basically determined by Y e and 
the entropies (and expansion time scales). In Fig. 25 we 
plot entropies (in ks per nucleon) and densities at the 
time of ejection for the three different neutron star spins 
(run E, I and J). It is very interesting to note that the 
densities at the moment of ejection are very high, the 
bulk of matter becomes unbound at densities from ~ 10 14 
to ~ 10 12 g cm~ 3 . This is well above the neutron drip 
(4 • 10 11 g cm -3 ), where in spite of the high temperatures 
very large (A ~ 250 according to the LS-EOS) and very 
neutron rich (Z/A ~ 0.15) nuclei are present in apprecia- 
ble amounts (Xh ~ 0.3). These nuclei are far from being 
experimentally well-known. Hence, to start r-process cal- 
culations from these initial conditions, very exotic nuclei 
(not in vacuo, but immersed in a dense neutron gas) have 
to be implemented in the corresponding reaction networks. 
In addition, the effects of the high Fermi-energies on the 
reaction rates including beta decays (Pauli-blocking) have 
to be accounted for. 

The temperatures are strongly dependent on the ejection 
mechanism which is closely related to the stellar spins. 
In the corotating runs the ejecta are initially located on 
the front side of each star (with respect to the orbital 
velocity). Due to gravitational torques this matter gets 
smoothly stretched and thereby decompressed, no spikes 
in pressure, temperature and the time dependent viscosity 
parameter a (see Eq. 12) are visible. In this case the tem- 
peratures are only slightly above the initial temperature 
(4 • 10 9 K; see Appendix A). We suspect this temperature 
increase to be mainly due to (artificially high) viscosity. 
The material is ejected in a different way for the other spin 
configurations. In the case without stellar spins the ejecta 
can be separated into two groups according to their ejec- 
tion mechanism. The material of the first group is found in 



the spiral arm structure, ejected in a way similar to corota- 
tion and is thus essentially cold. The second group comes 
from a region that gets strongly compressed and thereby 
heated up to temperatures around 6 MeV. In the follow- 
ing expansion, however, this material cools down quickly. 
In the case where the stars spin against their orbital mo- 
tion all the material is ejected by the second mechanism 
thereby reaching even higher temperatures in the com- 
pression phase (around 9 MeV) for a short time. 
Since the material that gets unbound in the coalescence 
of initially corotating systems stays essentially cold (10 8 
K, see Lai 1994), we expect the Y e of this matter to be 
close to the initial values of the cold neutron stars, i.e. 
0.01 < Y e < 0.05 with small contributions from the stellar 
crust. The cases with different stellar spins eject material 
that gets heated appreciably before being cooled by the 
expansion. In these cases Y e might be different from the 
initial values, since temperatures are high enough for the 
charged current reactions (e~— , e+— capture) to set in at 
non-negligible rates. 

Owing to the problems in explaining the observed r- 
process abundances entirely by type II supernovae, there 
seems to be a need for at least one further astrophysical 
scenario that is able to produce r-process nuclei in appre- 
ciable amounts. Neutron star mergers are attractive can- 
didates since they would in a natural way provide large 
neutron fluxes, low Y e s and moderate entropies (which 
provides r-process matter more easily than high Y e and 
entropy conditions). An r-process under such conditions 
should be very efficient and produce mostly elements in 
the high mass region. Thus, perhaps all of the r-process 
matter with A > 110, that can only be produced in the 
right amounts in supernova calculations if artificially high 
entropies are applied (see Freiburghaus et al. 1997; Taka- 
hashi et al. 1994), perhaps all of this matter could be 
synthesized in neutron star binary (or BH-NS) mergers. 
Assuming a core collapse supernova rate of 2.2- 10 -2 (year 
galaxy) - 1 (Ratnatunga 1989), one needs 10~ 6 to 10~ 4 M Q 
of ejected r-process material per supernova event to ex- 
plain the observations if type II supernovae are assumed 
to be the only source. The rate of neutron star mergers, 
which is by far more uncertain, has recently been esti- 
mated to be 8 • 10~ 6 (year galaxy) -1 (see van den Heuvel 
& Lorimer 1996). Taking these numbers, one would hence 
need ~ 3 • 10~ 3 M Q to <~ 0.3 M Q per event for an ex- 
planation of the observed r-process material exclusively 
by neutron star mergers. Thus our results for the ejected 
mass from 4 • 10~ 3 to 3 — 4 • 10~ 2 M Q look promising 
(see Fig. 26). Meyer (1989) calculated the decompression 
of initially cold neutron star material (T = K) assum- 
ing the expansion to be given by multiples of the free- 
fall time scale. He found that the decompressed neutron 
star material gives always, i.e. regardless of the expansion 
rate, rise to r-process conditions. Thus even the initially 
cold ejecta from corotating configurations should heat up 
during the expansion and form r-process nuclei. If, as sug- 
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gcsted by Meyer (1989), large parts of the ejected material 
should consist of r-process nuclei, neutron star mergers 
could account for the whole observed r-process material 
in the galaxy. However, whether the observed abundance 
patterns can be explained with this scenario remains an 
open question and is left to further investigations. 

4. Summary 

We have presented the results of three-dimensional New- 
tonian hydrodynamic simulations of the merger of equal 
mass neutron star binary systems. We followed the evolu- 
tion of two neutron stars containing 1.6 (1.4) M Q of bary- 
onic matter starting with an initial center of mass distance 
of 45 km for ~ 13 ms. In a total of 11 runs we tested the 
sensitivity of the results to the physical parameters of the 
binary system and to a number of model assumptions. 
In all cases we find a rapidly spinning central object (pe- 
riod ~ 1 ms) with masses (depending on the initial spins) 
between 2.5 and 3.1 M Q . This object might be stable on 
the simulation time scale, but we suspect it to collapse 
to a black hole within milliseconds after the merger. Our 
central objects are surrounded by thick disks containing 
(depending on the initial spins) ~ 0.1 to 0.3 M Q . In addi- 
tion we find long extended tails for the corotating models. 
These expand explosively when the neutron star matter 
gets decompressed by tidal torques in the LS-EOS cases 
when the adiabatic exponent rises and a few times 10 49 
erg are released due to the recombination of nucleons into 
heavy nuclei. For the polytropic EOS (r = 2.6) the tails 
remain thin and well-defined. For the other initial spins 
both the central object and the disk are embedded in a 
low density cloud of decompressed neutron star matter. In 
all cases (apart from the one where the stars spin against 
the orbit) almost baryon free funnels form above the poles 
of the central object. In these funnels large gradients of ra- 
diation pressure could be built up that would be able to 
accelerate indrifting matter into two jets pointing away 
from the poles. These funnels would also be an ideal place 
for relativistic fireballs to form. However, a mass as small 
as 10~ 5 Mq within such a fireball would be enough to 
prevent a GRB from forming. The present resolution is 
too low to draw conclusions on this point. We find SPH- 
smoothed temperatures of up to 50 MeV. These are found 
in macroscopic vortex structures that form along the con- 
tact surface of both stars and which we suspect to origi- 
nate from Kelvin- Helmholtz instabilities. 
The main new result is the amount of mass that is ejected 
into space. We find that, dependent on the initial spins, 
between 4 • 10~ 3 and 4 • 10~ 2 M become unbound for the 
realistic equation of state of Lattimer and Swesty. This 
result is strongly dependent on the EOS, a stiff polytrope 
ejects only around one half of this material. In the test 
case of a soft polytrope (r = 2.0) we cannot resolve any 
mass loss, which indicates a strong sensitivity on the stiff- 
ness of the EOS. The material gets ejected at very high 
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densities ranging from 10 12 to 10 14 gcm -3 . The bulk of 
matter gets ejected with Ye below 0.05 with small con- 
taminations of the neutron star crust (Y e w 0.3). Such 
low Y~ e , low entropy matter is prone to form r-process nu- 
clei. However, the results concerning Y e are biased by our 
simple neutrino treatment. Using recent rates for neutron 
star mergers we find that ~ 3 • 10~ 3 to 0.3 M Q of r-process 
material would have to be ejected to explain the observed 
abundances exclusively by coalescing neutron stars. Thus 
our numbers for the amount of ejecta look promising and 
if, as suggested, large parts of this matter consist of r- 
process nuclei, neutron star mergers could account for all 
the observed r-process material in the Galaxy. 
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A. The Particle Setup 

In the initial setup of the SPH-particles we adopt an inter- 
mediate approach between the extremes of total randomness 
and the well defined order of a lattice. In a first step we create 
a large set of spherical shells, containing initially randomly 
distributed particles which have been relaxed according to 
some Coulomb-type force law. The so constructed shells are 
then joined in a concentric way, and the particles are assigned 
masses that, together with the smoothing length, reproduce 
the desired density profile. This configuration is finally relaxed 
with the SPH-code. 

To construct a single shell of a given particle number we dis- 
tribute the particles randomly on the surface of a unit sphere. 
We assume some repulsive, dissipative Coulomb-type force law: 

/.=E§-^ ( A1 ) 

where a is some dissipation parameter, r;j is the difference 
vector of particle i and j, rij = t% — rj, r t and Vi are particle 
position and velocity vectors. In the following relaxation pro- 
cess only the force component tangential to the spherical sur- 
face was used in order to constrain the movement to a constant 
radius. After each integration time step the particle positions 
were projected back onto the surface. 

To be able to vary the particle number density with the radius 
we introduce a parameter < e < 1. The distance to the new 
shell Ar n+i is then calculated from the previous one Ar n via 

Ar n+1 = eAr„. (A2) 

Denoting the distance from the center to the first shell by 
An = S and requiring that 

iJstar = Ar " ( A3 ) 
i=l 
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-Rstar is the radius of our star, we find for a given number of 
shells -/Vsheii and given e ^ 1: 



8 = R B 



(A4) 



Now the question arises of how to choose the particle number 
in a given shell. We have as a constraint that typical particle 
separations within shell n should be the same as typical dis- 
tances between shells Af n = Ar " + ^ r " +1 . We consider now a 
sphere of radius R with iV n SPH-particles distributed on its 
surface. If we calculate the typical distance of two particles 
within one shell and equate it to Af„ we find for the particle 
number in shell n 



N„ 



4.R 2 



In the case of N n > 1 we have 
expand the square root in the denominator: 



R 



(A5) 



<C 1 and can therefore 



N n « 16 



R 2 



(A6) 



The same result is recovered by assigning to each particle the 
surface of a circle and setting this equal to . 
To find the optimal relative orientation of the shells we go back 
to the idea of the Coulomb-type forces. 

The innermost shell is kept fixed. Then the torque exerted 
from shell one to shell two is calculated and shell two is ro- 
tated (according to the corresponding Eulerian angles) until 
the torque vanishes. Fixing the second shell at the new posi- 
tion, the torque on shell three from the inner shells is calculated 
and so on. 

The profiles p p (r), Y p (r) and u p (r) define our initial neutron 
star model. They are calculated using a Newtonian ID stel- 
lar structure code and a self-consistent table of p, Y e and u 
generated using the LS-EOS. For the calculation of Y e , we use 
(C6) (at T — 4 • 10 9 K) and the [3- equilibrium condition for the 
chemical potentials: 



Me = A + Q — Pi> e 



(A7) 



where Q is the mass difference of neutron and proton times c 
and the difference in the chemical potentials of neutron and 
proton (without rest masses), fi, is provided by the LS-EOS. 
Since we want to construct a cold neutron star where no neutri- 
nos are present, we assume their chemical potentials to vanish, 



fJ,p e = 0. 



(A8) 



The y e -profiles are for numerical reasons restricted to values 
^ 0.05. 

After values of Y e and u have been assigned to the particles 
according to the profiles, the masses still have to be distributed. 
In principle 



p — Wm, 



(A9) 



has to be solved, where the components of the vectors are the 
values at the particle positions (density and mass) and W is 
the kernel matrix. 

In order to enforce sphericity, we assign one single value for 



all the masses mi and one for all the smoothing lengths hi 
to all the particles in a given shell. We prescribe an allowed 
range of neighbours for the mean neighbour number in each 
shell (the neighbours are found by traversing the tree of our 
SPH-code) and adjust the smoothing lengths correspondingly 
(typical neighbour numbers are ?a 65). Starting with the masses 
from the density self-contribution of each particle 



nn = irhi p p (ri), 



(MO) 



where the spherical spline kernel (see Benz (1990)) has been 
used, the masses are iterated until 



max[ \£i^Ml) <ep , 

j V P p (rj) ! 



(All) 



where pj is the mean value of the density in shell j (e p was 
typically a few times 10 -3 ). 

This particle setup was then finally relaxed with the SPH-code 
where a velocity dependent additional force term was applied: 



fi, hydro > fi, hydro fVi. 



(A12) 



The parameter 7 gives a characteristic damping time scale 
t = 7 _1 . To obtain an efficient convergence towards the equi- 
librium solution it is important to choose 7 to be slightly over- 
critical, i.e. t < Tosc, where the typical oscillation period Tosc IS 
approximately given by the sound crossing time of the neutron 
star r osc ~ t s — (Gp) _ 2 0.3 ms. Some properties of the 
relaxed star are shown in Fig. 27. 

B. The Equation of State 

For all our calculations the default set of nuclear parameters 
(see Lattimer & Swesty (1991)) with K = 180 MeV is used. 
Most of the scheme that we are going to describe here has 
been developed for version 2.6 of the LS-EOS. We encountered 
several regions (especially for low Y e and densities where large 
amounts of the heavy nuclei are present) where the Newton- 
Raphson iteration to solve the set of equilibrium equations did 
not converge for the default set of guess values (nucleon number 
density inside nuclei m , nucleon degeneracy parameters r\ n and 
r/p) due to very small convergence radii. 

To find accurate guess values gi in the above mentioned critical 
regions we use the fact that the gi are slowly varying functions 
of Y e . We chose a fixed, uncritical value of Y e , Y^ c{ = 0.25, for 
which reference surfaces of the guess values G i Y ret over the n — 
T-plane, n is the nucleon number density, T the temperature, 
were calculated 



G.yref : (n,T)-» 9i . 



(Bl) 



To patch the surface over the critical region for some given Y e , 
the corresponding rectangular piece in the reference surface 
was cut out and matched to the hole in the surface of the 
desired Y e . Let 



(ni,n 2 ) x (Ti,T 2 ) 



(B2) 



be the region where the code does not converge for a given Y e . 
We then calculate multipliers rrn,j for each guess value gi 



9i,j(Ye) 



(B3) 
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with j labeling the edge point, in such a way that the edge 
points in the reference surface coincide with those for the crit- 
ical Y e when multiplied with the corresponding rriij: 



where a « 1.25, E Ve is the energy of the emitted neutrino, 



9i,j(Y e ) = mi,j(Y e ) gi,j(Y* e ). 



(B4) 



Then a multiplier rm for the points within the critical region 
is calculated by means of linear interpolation: 



m t (n,T,Y e ) = (1 - x)(l - y) mi,i(Y e ) + x (1 - y) m t , 2 (Y e ) 
+ xy mi, 3 (Y e ) + (1 - x)y m iA (Y e ), (B5) 

where x = (n - m)(n 2 - m)" 1 and y = (T - Ti)(T 2 - Ti)" 1 
and the edge points are numbered starting with (m,Ti) in a 
counter-clockwise sense. The guess value variables in the criti- 
cal region are then approximated by 



9i(n, T, Y e ) = mi(n, T, Y e ) ■ g^n, T, Y e rc ). 



(B6) 



Using this scheme accurate tables for the guess values have 
been calculated. Using these tables the code converges rapidly 
and safely everywhere. 

In our hydrodynamic calculations we use a tabular form of the 
LS-EOS (version 2.7), where the above described scheme has 
been used. Our table contains 153 entries in log(p), 121 entries 
in log(w), and 25 entries in Y e . For our calculations we need 
to tabulate the pressure, the temperature and the difference 
in the nucleon chemical potentials, fi. The abundances are dif- 
ficult to tabulate since they may vary enormously from grid 
point to grid point. Therefore, every time that abundances are 
needed (e.g. to calculate the total nuclear binding energy of 
the system), the original LS-EOS is used. 

Since we use the specific internal energy as independent vari- 
able, we performed a bisection iteration to find the correspond- 
ing temperatures. If non-monotonic values in T were encoun- 
tered, the iteration used u values that were averaged over 
neighboring values. Note, however, that apart from that, no 
other manipulation of the EOS-data has been performed (Ruf- 
fert et al. (1996) used smoothed EOS-tables). 
The interpolation in the table is a delicate topic in itself since 
thermodynamic consistency for the values between the tabu- 
lated grid points is not guaranteed. An inconsistent interpola- 
tion, i.e. an interpolation that does not fulfill the constraints 
posed by the Maxwell relations of thermodynamics, will lead to 
an artificial buildup of entropy, or temperature. Swesty (1996) 
has proposed an interpolation scheme insuring thermodynamic 
consistency as well as the continuity of the derivatives of the 
thermodynamic functions. However, this approach applied to 
our 3D-input parameter space would require 216 terms and 
the evaluation of fourth order derivatives for each table call. 
Clearly, this procedure is - at least at present - computation- 
ally prohibitive. Considering the other approximations in our 
model, we think that we can justify just a linear interpolation 
in our table. 

C. Neutrino Emission Rates 

Starting from the electron capture cross section (see Tubbs and 
Schramm (1975)), which reads with our approximations 



<?EC 



1 + 3^ 



O~0 



(m e c 2 ) 2 



(CI) 



<7<) = 



(C2) 



and Gf the Fermi constant, the number of electron captures 
per time and volume, Rec, is given by (k B = 1) 



Rec = Vpn 



* !+3Q aoT 5 ^]. 



h s c 2 (m e c 2 ) 2 



(C3) 



The corresponding energy loss in neutrinos per time and vol- 
ume is 



Qec = ri vn 



7T 1 + 3 a 2 
h 3 c 2 (m e c 2 ) 2 



o T 6 F 5 [ Ve 



(C4) 



Here rj pn = pNA(Y n — Y p )/(e r — l) ; with the Yi denoting 
neutron and proton abundances, [n the corresponding chemical 
potentials (without rest mass) and Na Avogadro's constant. 
The factor rj pn , derived by Bruenn (1985) under the assumption 
that 4-momentum transfer between leptons and nucleons is 
negligible ( "elastic approximation" ) , takes into account effects 
resulting from the phase space restrictions of the final state 
nucleons. F5 is the usual Fermi integral given by 



F„(z) 



f 

Jo 



+ 1 



-da;. 



(C5) 



rj e is the electron degeneracy parameter, ij e = ^f-- The electron 
chemical potential is calculated according to 

fi e = kT(pN A )^ {[\/b 2 +a3 + b]i - Wb 2 + a? - 6]*}, (C6) 

where a = (tt 2 /3 - P 2 /2){pN A )- 2 '* , b = f^(#) 3 y e , /J = 
m e c 2 T 1 have been introduced (see Baron (1985)). 
Assuming p e + — —p e the corresponding formulae for the 
positron captures read 



Rpc — r\ nv 



and 



7T 1 + 3 a 2 
hFc* (rn e c 2 ) 2 



ao T 5 Fi[-i]e] 



(C7) 



(C8) 



where rj np is found from r) pn by interchanging the indices. For 
simplicity we assume here that matter consists only of protons 
and neutrons: 



Yp — Y e , Y n — 1 Y e . 



The rate of change in the electron fraction is given by 
• _ Rpc y _ Rec y 



(C9) 



(CIO) 



16 



S. Rosswog et al.: Mass ejection in neutron star mergers 



References 

Abramovici, A., et al. 1992, Science, 256, 325 
Bardeen, J.M., Press, W.H., & Teukolsky, S.A. 1972, ApJ, 178, 
347 

Baron, E. 1985, Ph.D. thesis, State University of New York at 
Stony Brook 

Baumgarte, T.W., Cook, G.B., Scheel, M.A., Shapiro, S.L., & 
Teukolsky, S.A. 1997, Phys. Rev. Lett., 79, 1182 

Benz, W. 1990, in Numerical Modeling of Stellar Pulsations, 
ed. J. Buchler (Dordrecht: Kluwer Academic Publishers), 
269 

Benz, W., Bowers, R.L., Cameron, A.G.W., & Press, W.H. 

1990, ApJ, 348, 647 
Bildsten, L., & Cutler, C. 1992, ApJ, 400, 175 
Bradaschia, C, et al. 1990, Nucl.Instrum. Methods Phys. Res. 

A, 289, 518 
Bruenn, S. W. 1985, ApJS, 58, 771 
Costa, E., et al. 1997, Nature, 387, 783 

Davies, M.B., Benz, W., Piran, T., & Thielemann, F.-K. 1994, 
ApJ, 431, 742 

Djorgovski, S.G., et al. 1997, Nature, 387, 876 

Drazin, P., & Reid, W. 1981, Hydrodynamic Stability (Cam- 
bridge: Cambridge University Press) 

Eichler, D., Livio, M., Piran, T., & Schramm, D. N. 1989, Na- 
ture, 340, 126 

Flowers, E., & Itoh, N. 1979, ApJ, 230, 847 

Frail, D.A., et al. 1997, Nature, 389, 261 

Freiburghaus, C, Rauscher, T., Thielemann, F.-K., Kratz, K.- 
L., & Pfeiffer, B. 1997, in Proc. 4th Int. Symp. on Nu- 
clei in the Cosmos (Univ. of Notre Dame), ed. J. Gorres, 
G. Mathews, S. Shore, & M. Wiescher, Vol. A621, Nuc. 
Phys., 405c 

Freiburghaus, C, Rembges, J.-F., Rauscher, T., Kolbc, E., 
Thielemann, F.-K., Kratz, K.-L., Pfeiffer, B. & Cowan, J.J. 
1998, ApJ submitted 

Galama, T., et al. 1997, Nature, 387, 479 

Goodman, J., 1986, ApJ, 308, L51 

Haft, M., Raffelt, G., & Weiss, A. 1994, ApJ, 425, 222 

Hernquist, L., & Katz, N. 1989, ApJS, 70, 419 

Itoh, N., et al. 1989, ApJ, 339, 354 

Klebesadel, R.W., Strong, I.B., & Olsen, R.A. 1973, ApJ, 182, 
L85 

Kochanek, C.S. 1992, ApJ, 398, 234 

Lai, D. 1994, MNRAS, 270, 611 

Lai, D. 1996, Phys.Rev.Lett., 76, 4878 

Lai, D., Rasio, F.A., & Shapiro, S.L. 1993a, ApJS, 88, 205 

Lai, D., Rasio, F.A., & Shapiro, S.L. 1993b, ApJ, 406, L63 

Lai, D., Rasio, F.A., & Shapiro, S.L. 1994a, ApJ, 420, 811 

Lai, D., Rasio, F.A., & Shapiro, S.L. 1994b, ApJ, 423, 344 

Lai, D., Rasio, F.A., & Shapiro, S.L. 1994c, ApJ, 437, 742 

Lai, D., & Shapiro, S.L. 1995, ApJ, 443, 705 

Landau, L., & Lifshitz, E. 1976, Mechanik (12. ed.) (Berlin: 

Akademieverlag) 
Lattimer, J. M., & Schramm, D. N. 1974, ApJ, (Letters), 192, 

L145 

Lattimer, J. M., & Schramm, D. N. 1976, ApJ, 210, 549 
Lattimer, J. M., & Swesty, D. 1991, Nucl. Phys., A535, 331 
Lattimer, J. M., & Yahil, A. 1989, ApJ, 340, 426 
Lipunov, V. et al. 1995, ApJ, 454, 593 
Lombardi, J., & Rasio, F.A. 1997, Phys. Rev., D56, 3416 



Luck, H. 1997, Class. Quant. Grav., 14, 1471 
Mathews, G.J., & Wilson, J. R. 1997, ApJ, 482, 929 
Metzger, M.R., et al. 1997, Nature, 387, 878 
Meyer, B. S. 1989, ApJ, 343, 254 

Meyer, B. S., & Brown, J. S. 1997, in Proc. 4th Int. Symp. on 
Nuclei in the Cosmos (Univ. of Notre Dame), ed. J. Gorres, 
G. Mathews, S. Shore, & M. Wiescher, Vol. A621, Nuc. 
Phys., 409c 

Misner, C, Thorne, K., & Wheeler, J. 1973, Gravitation (New 

York: Freeman) 
Monaghan, J.J. 1992, Ann. Rev. Astron. Astrophys., 30, 543 
Monaghan, J.J., & Gingold, R. 1983, J. Comp. Phys., 52, 374 
Morris, J.P., & Monaghan, J.J. 1997, J. Comp. Phys., 136, 41 
Miiller, E. 1998, in Computational Methods for Astrophysical 

Fluid Flow, ed. O. Steiner & A. Gautschy (Berlin: Springer 

Verlag) 

Narayan, R., Piran, T., & Shemi, A. 1991, ApJ, 379, L17 

Narayan, R., Paczyhski, B., & Piran, T. 1992, ApJ, 395, L83 

Paczyhski, B. 1986, ApJ, 308, L43 

Paczyhski, B. 1990, ApJ, 363, 218 

Paczyhski, B. 1991, Acta Astron., 41, 257 

Paczyhski, B. 1992, Nature, 355, 521 

Peters, P., & Mathews, J. 1963, Phys. Review, 131, 435 

Peters, P., & Mathews, J. 1964, Phys. Review, B136, 1224 

Piran, T., & Shemi, A. 1993, ApJ 403, L67 

Prakash, M. 1997, in Nuclear and Particle Astrophysics, ed. 

J. Hirsch & D. Page (Cambridge: Cambridge University 

Press) 

Prakash, M., Cooke, J.R., & Lattimer, J.M. 1995, Phys. Rev., 
D52, 661 

Rasio, F.A., & Shapiro, S.L. 1992, ApJ, 401, 226 
Rasio, F.A., & Shapiro, S.L. 1994, ApJ, 432, 242 
Rasio, F.A., & Shapiro, S.L. 1995, ApJ, 438, 887 
Ratnatunga, K.U., & van den Berg, S. 1989, ApJ, 343, 713 
Rees, M., & Mesczaros, P. 1992, MNRAS, 258, 41 
Rosswog, S., Liebendorfer, M., Thielemann, F.-K., Davies, 
M.B., Benz, W., & Piran, T. 1998a, in Proceedings of the 
2 nd Oak Ridge Symposium on Atomic and Nuclear Astro- 
physics 

Rosswog, S., Liebendorfer, M., Thielemann, F.-K., Davies, 
M.B., Benz, W., & Piran, T. 1998b, in Proc. 9th Work- 
shop on Nuclear Astrophysics, Schloss Ringberg, Germany, 
March 1998 

Ruffert, M., Janka, H.-T., & Schafer, G. 1996, A & A, 311, 532 
Ruffert, M., Janka, H.-T., Takahashi, K., & Schafer, G. 1997, 

A & A, 319, 122 
Ruffert, M., Rampp, M., & Janka, H.-T. 1997, A & A, 321, 991 
Sahu, K.C., et al. 1997, Nature, 387, 476 
Schinder, P.J., et al. 1987, ApJ, 313, 531 

Shapiro, S.L., & Teukolsky, S. A. 1983, Black Holes, White 
Dwarfs and Neutron Stars (New York: Wiley & Sons) 

Shemi, A., & Piran, T. 1990, ApJ 379, L17 

Shibata, M., Nakamura, T., & Oohara, K. 1992, Prog. Theor. 
Phys. 88, 1079 

Shibata, M., Nakamura, T., & Oohara, K. 1993, Prog. Theor. 
Phys. 89, 809 

Shibata, M., Baumgarte, T. W., & Shapiro, S.L. 1998, preprint, 

gr-qc/9805026 
Stark, R.F., Piran, T. 1985, Phys. Rev. Lett., 55, 891 
Swesty, F. 1996, J. Comput. Phys., 127, 118 



S. Rosswog et al.: Mass ejection in neutron star mergers 



Symbalisty, E. M. D., & Schramm, D. N. 1982, Astrophys. 
Lett., 22, 143 

Takahashi, K., Witti, J., & Janka, H.-T. 1994, A & A, 286, 857 
Taylor, J.H. 1994, Rev. Mod. Phys., 66, 711 
Thorne, K. S. 1997, preprint, gr-qc/9706057 
Thorsett, S.E. 1996, Phys. Rev. Lett., 77, 1432 
Tubbs, D., & Schramm, D. 1975, ApJ, 201, 467 
van den Heuvel, E. 1991, in Neutron Stars: Theory and Ob- 
servation, eds. Ventura, J. & Pines, D. (Dordrecht: Kluwer 
Academic Publishers) 
van den Heuvel, E. 1994, in Interacting Binaries, eds. Nuss- 

baumer, H. & Orr, A. (Berlin: Springer Verlag) 
van den Heuvel, E., & Lorimer, D. 1996, MNRAS, 283, L37 
van Paradijs, J. V., et al. 1997, Nature, 386, 686 
Weber, F., & Glendenning, N.K. 1992, ApJ, 390, 541 
Weber, F., Glendenning, N.K., & Weigel, M.K. 1991, ApJ, 373, 
579 

Weber, F., & Weigel, M.K. 1989, Nuc. Phys., A505, 779 
Wilson, J. R., & Mathews, G.J. 1995, Phys. Rev. Lett., 75, 
4161 

Wilson, J. R., Mathews, G.J., & Marronetti, P. 1996, Phys. 

Rev., D54, 1317 
Wiseman, A.G. 1997, Phys. Rev. Lett., 79, 1189 
Zhuge, X., Centrella, J., & McMillan, S. 1994, Phys. Rev., D50, 

6247 

Zhuge, X., Centrella, J., & McMillan, S. 1996, Phys. Rev., D54, 
7261 



18 



S. Rosswog et al.: Mass ejection in neutron star mergers 





Fig. 1. Morphology of run A (representative for corotation). 
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Fig. 2. Morphology of run C (corotation, polytropic EOS with 
T = 2.6); the large dots denote the positions of the escaping 
particles. 
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Fig. 3. Morphology of run I (no initial spin). 
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Fig. 4. Morphology of run J (spins against the orbital motion). 




Fig. 5. Total amount of nuclear energy (in units of 10 50 erg) 
in the system (time is given in ms). 
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Fig. 6. The panels show the evolution of total angular mo- 
mentum L and total energy E of runs A, I and J. L and E 
are conserved to approximately 3 ■ 10~ 3 (the decrease in the 
beginning corresponds to the phase where angular momentum 
as well as energy are lost through the emission of gravitational 
waves) . 
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Fig. 7. Cut through the y-z-plane of the last dumps of the 
runs E (corotation, no neutrinos) and G (corotation, neutrinos 
in free-streaming limit); the labels in all density contour plots 
refer to log(p [gcm - 3] ). 
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Fig. 8. Last dumps of the runs I (no initial spins) and J (spins 
against orbit). 
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Fig. 9. The figure shows the mass that has a larger density 
than log(p) for all the runs (last dump). 
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Fig. 10. Distribution of mass with cylindrical radius (last 
dump). Masses are given in solar units, the radius in km. 
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Fig. 11. Time evolution (run A, corotation) of the mass that 
has a larger density than log(p). The contour lines in the 
t-log(p)-plane correspond to 2.8, 2.9, 3.0 and 3.1 Mq. Time 
is given in ms, p in g cm" 3 . 
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Fig. 12. Density contours of run E (corotation), lengths are 
given in km. 
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Fig. 13. Density contours of run I (no initial spins), lengths 
are given in km. 
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Fig. 14. Density contours of the runs E (LS-EOS, left column) 
and C (polytrope with F = 2.6, right column), lengths are given 
in km, the contour labels refer to log(p[ gcm -3]). 
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Fig. 15. Cut through the x-y-plane of the runs I (no spins, left 
column) and J (spins against orbit, right column), lengths are 
given in km. 
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Fig. 16. Maximum densities (in units of gem 3 ) obtained in 
our runs. 
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Fig. 17. The upper three curves correspond to the Ke- 
pler-velocities (c = 1), the lower ones are mean tangential 
particle velocities versus cylindrical radius. The solid line cor- 
responds to corotation (run E), the circles to run I (no spins) 
and the asterisks to run J (spins against orbit). 
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Fig. 18. Shown are the maximum temperatures of the different 
runs. The upper panel shows the maximum particle tempera- 
tures, the lower one the maximum of the SPH-smoothed val- 
ues. The legend refers to both panels, temperatures are given 
in units of MeV. 
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Fig. 19. The left column shows the SPH-smoothed tempera- 
tures (in MeV) of run E (corotation). The right column shows 
the origin of all particles from either of the stars (crosses, star 
1; dots, star 2). 
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Fig. 20. SPH-smoothed temperatures (in MeV) of run I (no 
initial spins). The particles from each star are plotted with 
different symbols. 
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Fig. 21. Same as Fig. 20, but for run J (spins against orbital 
angular momentum). 
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Fig. 22. Comparison of nuclear binding energy and the total 
energy emitted in neutrinos for the three morphological regions 
of the coalesced object. The upper panel refers to the central 
object, the one in the middle to the disk and the lowest to the 
spiral arms. On the abscissa time (in ms) and on the ordinate 
the logarithm of the energies (in erg; B nuc and E„, t ot denote nu- 
clear binding energy and the total emitted neutrino energy) is 
shown. In the central object (in this simple model) the neutrino 
emission dominates clearly over the released nuclear binding 
energy. The "bump" in panel one short before contact results 
from nuclei that form when the density decreases due to tidal 
stretching. In the low mass spiral arms (see Table 2) clearly 
the nuclear energy deposition (~ 8 • 10 49 erg) dominates. 




Fig. 23. Mass ejection (in solar units) of the different models. 
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Fig. 24. Logarithm of total nuclear binding energy (in erg) 
present in the escaping mass. 
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Fig. 25. Entropies and densities of the particles at the moment 
of ejection. Dots refer to run E (corotation), open circles to run 
I (no spin) and asterisks to run J (spins against orbit). 
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Fig. 26. The shaded region shows the amount of ejected mate- 
rial found in our calculations. The circle shows the amount of 
ejecta per event if SN II are assumed to be the only sources of 
the r-process. The asterisk gives the needed ejecta per merging 
event for the rate of Narayan et al. (1991), the cross for the 
estimate of van den Heuvel and Lorimer (1996). The event rate 
is given in year -1 galaxy -1 , the ejected mass in solar units. 
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Fig. 27. Properties of the configuration after relaxation. The 
upper left (upper right, lower right) panel shows the distribu- 
tion of the particle masses (densities, smoothing lengths) with 
radius, the lower left panel shows the neighbour numbers of 
each particle. 
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